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We present a comprehensive theoretical study of the magnetic held dependence of the near-held 
radiative heat transfer (NFRHT) between two parallel plates. We show that when the plates are 
made of doped semiconductors, the near-held thermal radiation can be severely affected by the 
application of a static magnetic held. We hnd that irrespective of its direction, the presence of a 
magnetic held reduces the radiative heat conductance, and dramatic reductions up to 700% can be 
found with helds of about 6 T at room temperature. We show that this striking behavior is due 
to the fact that the magnetic held radically changes the nature of the NFRHT. The held not only 
affects the electromagnetic surface waves (both plasmons and phonon polaritons) that normally 
dominate the near-held radiation in doped semiconductors, but it also induces hyperbolic modes 
that progressively dominate the heat transfer as the held increases. In particular, we show that 
when the held is perpendicular to the plates, the semiconductors become ideal hyperbolic near-held 
emitters. More importantly, by changing the magnetic held, the system can be continuously tuned 
from a situation where the surface waves dominate the heat transfer to a situation where hyperbolic 
modes completely govern the near-held thermal radiation. We show that this high tunability can be 
achieved with accessible magnetic helds and very common materials like n-doped InSb or Si. Our 
study paves the way for an active control of NFRHT and it opens the possibility to study unique 
hyperbolic thermal emitters without the need to resort to complicated metamaterials. 


I. INTRODUCTION 

Radiative heat transfer is a topic with numerous funda¬ 
mental and technological implications across disciplines 
[T]. Presently, the investigation of radiative heat trans¬ 
fer between closely spaced objects is receiving a great 
attention EHS]- The basic challenges now are the under¬ 
standing of the mechanisms governing thermal radiation 
at the nanoscale and the ability to control this radiation 
for its use in novel applications. During a long time, it 
was believed that the maximum radiative heat transfer 
between two objects was set by the Stefan-Boltzmann 
law for black bodies. However, this is only true when ob¬ 
jects are separated by distances larger than the thermal 
wavelength (9.6 /im at room temperature) and the ra¬ 
diative heat transfer is dominated by propagating waves. 
When two objects are brought in closer proximity, the 
thermal radiation is dominated by interference effects 
and, more importantly, by the near field emerging from 
the materials surfaces in the form of evanescent waves. 
This was first established theoretically by Polder and 
Van Hove in 1971 [ 7 ] using Rytov’s framework of fiuc- 
tuational electrodynamics These researchers pre¬ 

dicted that the NFRHT could overcome the black body 
limit by several orders of magnitude, achieving what is re¬ 
ferred to as super-Planckian thermal emission. Although 
this NFRHT enhancement was already hinted in several 
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experiments in the late 1960 ’s [iniiii], its unambiguous 
confirmation came only in recent years [UHll]. These 
experiments, in turn, have triggered off the hope that 
NFRHT may have an impact in different technologies 
such as heat-assisted magnetic recording [28l [29] , ther¬ 
mal lithography m, scanning thermal microscopy isi- 
l33] . coherent thermal sources [3ll |35], near-field based 
thermal management [22l l3Hl438j . thermophotovoltaics 
ISIIIO] and other energy conversion devices m- 

Presently, one of the central research lines in the field 
of radiative heat transfer is the search for materials where 
the NFRHT enhancement can be further increased. So 
far, the largest enhancements have been experimentally 
reported in polar dielectrics [igiTniiii], where the near¬ 
field thermal radiation is dominated by the excitation 
of surface phonon polaritons (SPhPs) |42]. Similar en¬ 
hancements have been predicted and observed in doped 
semiconductors due to surface plasmon polaritons (SPPs) 
[23l |27l |43l |44|. Recently, it has been predicted that 
hyperbolic metamaterials could behave as broadband 
super-Planckian thermal emitters [45HT7] . Hyperbolic 
metamaterials are a special class of highly anisotropic 
media that have hyperbolic dispersion. In particular, 
they are uniaxial materials for which one of the principal 
components of either the permittivity or the permeability 
tensors is opposite in sign to the other two principal com¬ 
ponents [48]. Hyperbolic media have been mainly real¬ 
ized by means of hybrid metal-dielectric superlattices and 
metallic nanowires embedded in a dielectric host gano]. 
It has been demonstrated that these metamaterials ex¬ 
hibit exotic optical properties such as negative refrac- 
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tion, subwavelength imaging and focusing, and they can 
be used to do density of states engineering milO]. In 
the context of radiative heat transfer, what makes these 
metamaterials so special is the fact that they can support 
electromagnetic modes that are evanescent in a vacuum 
gap, but they are propagating inside the material. This 
leads to a broadband enhancement of transmission effi¬ 
ciency of the evanescent modes [46]. This special prop¬ 
erty has motivated a lot of theoretical work on the use 
of hyperbolic metamaterials for NFRHT mHH]. How¬ 
ever, no experimental investigation of this issue has been 
reported so far, which is mainly due to the difficulties 
in handling these metamaterials. In this sense, it would 
be highly desirable to find much simpler realizations of 
hyperbolic thermal emitters and, ideally, with tunable 
properties. 

Another key issue in the field of radiative heat transfer 
is the active control and modulation of NFRHT. In this 
respect, several proposals have been put forward recently. 
One of the them is based on the use of phase-change ma¬ 
terials (62] |63| , where the change of phase leads to a sig¬ 
nificant change in the material dielectric function. These 
materials include an alloy called AIST, where the phase 
change can be induced by applying an electric field [63] . 
and VO 2 , which undergoes a metal-insulator transition 
as a function of temperature [62|. It has also been sug¬ 
gested that the NFRHT between chiral materials with 
magnetoelectric coupling can be tuned by ultrafast op¬ 
tical pulses M- Another proposal to tune the NFRHT 
is to use ferroelectric materials under an external electric 
field [65| , although the predicted changes are rather mod¬ 
est (< 17%). Let us also mention that very recently it 
has been proposed that the heat flux between two semi¬ 
conductors can be controlled by regulating the chemical 
potential of photons by means of an external bias [66] . So 
in short, although these proposals are certainly interest¬ 
ing, some of them are not easy to implement and others 
are either not very efficient or they are restricted to very 
specific materials. In this sense, the challenge remains 
to introduce strategies to actively control NFRHT in an 
easy and relatively universal way. 

In this work we tackle and resolve some of the open 
problems described above by presenting an extensive the¬ 
oretical analysis of the influence of an external dc mag¬ 
netic field in the radiative heat transfer between two par¬ 
allel plates made of a variety of materials. We show that 
an applied magnetic field can indeed largely affect the 
NFRHT in a broad class of materials, namely doped (po¬ 
lar and non-polar) semiconductors. We find that, irre¬ 
spective of its orientation, the magnetic field reduces the 
NFRHT with respect to the zero-field case and we show 
that the reduction can be as large as 700% for fields of 
about 6 T at room temperature. This effect originates 
from the fact that the magnetic field not only strongly 
modifies the surface waves that dominate the NFRHT 
in doped semiconductors (both SPhPs and SPPs), but 
it also generates broadband hyperbolic modes that tend 
to govern the heat transfer as the field is increased. In 


particular, when the applied field is perpendicular to 
the plates surfaces, the semiconductors behave as hyper¬ 
bolic thermal emitters with highly tunable properties. By 
changing the field magnitude one can continuously tune 
the system and realize situations where (i) surface waves 
dominate the NFRHT, (ii) both surface waves and hyper¬ 
bolic modes contribute significantly to the near-held ther¬ 
mal radiation, and (iii) only hyperbolic modes contribute 
to the NFRHT and surface waves cease to exist. On the 
other hand, when the held is parallel to the surfaces the 
NFRHT is nonmonotonic as a function of the magnetic 
held. For moderate helds, surface waves and hyperbolic 
modes coexisting, while for high helds the NFRHT is 
largely dominated by hyperbolic modes. We emphasize 
that all these striking predictions are amenable to mea¬ 
surements and do not require the use of any complicated 
metamaterial. Thus, our work ohers a simple strategy to 
actively control NFRHT in a broad variety of materials 
and it also provides a very appealing recipe to realize hy¬ 
perbolic materials and, in particular, hyperbolic thermal 
emitters with highly tunable properties. 

The remainder of this paper is structured as follows. 
Section [H] describes the system under study and the gen¬ 
eral formalism for the description of NFRHT in the pres¬ 
ence of a magnetic held. We then turn in Sec. m to the 
application of the general results to the case of n-doped 
InSb as an example of a polar semiconductor. We discuss 
in this section both the results for diherent magnetic held 
orientations and the realization of highly tunable hyper¬ 
bolic thermal emitters. Section UTII is devoted to the case 
of Si as an example of non-polar semiconductor. Section 
IV summarizes our main results and discusses future di¬ 
rections. Finally, four appendixes contain the technical 
details of the general formalism and some additional cal¬ 
culations that support the claims in this paper. 


II. RADIATIVE HEAT TRANSFER IN THE 
PRESENCE OF A MAGNETIC FIELD: 

GENERAL FORMALISM 

Our main goal is to compute the radiative heat trans¬ 
fer in the presence of an external dc magnetic held within 
the framework of huctuational electrodynamics 0 ig. 
For simplicity, we shall concentrate here in the heat ex¬ 
changed between two inhnite parallel plates made of ar¬ 
bitrary non-magnetic materials and that are separated 
by a vacuum gap of width d, see Fig. The magnetic 
held can point in any direction and following Fig. [^a), 
we shall refer to the left plate as medium 1, the vacuum 
gap as medium 2, and the right plate as medium 3. 

When a magnetic held is applied to any object, it re¬ 
sults in an optical anisotropy that can be described by 
the following general permittivity tensor m 
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FIG. 1. (a) Schematic representation of the system under 

study: two parallel plates at temperatures Ti and T 3 sepa¬ 
rated by a vacuum gap of width d. (b) Heat transfer coeffi¬ 
cient of n-doped InSb as a function of the gap at zero mag¬ 
netic field. We show the total result and the individual con¬ 
tributions of s- and p-polarized waves (both propagating and 
evanescent), (c) The corresponding zero-field spectral heat 
flux as a function of frequency (and wavelength) for three 
different gaps. 


where according to Fig. [^a), x and y lie in the interface 
planes and ^ corresponds to the surface normal. The 
components of the permittivity tensor depend on the ap¬ 
plied magnetic field, as we shall specify below, and on the 
frequency (local approximation). Let us recall that the 
off-diagonal elements in Eq. Q are responsible for all 
the well-known magneto-optical effects (Faraday effect, 
Kerr effects, etc.) [67]. Thus, our problem is to compute 
the radiative heat transfer between two anisotropic par¬ 
allel plates. This generic problem has been addressed by 


Biehs et al [68] and we just recall here the central re¬ 
sult. The net power per unit of area exchanged between 
the parallel plates is given by the following Landauer-like 
expression [68] 

Q = ^°°^[0iH-©3M] J J^riu;,k,d), (2) 

where 0i(co’) = huj/[ex.-p{huj/kBTi) — 1], Ti is the absolute 
temperature of the layer i, cj is the radiation frequency, 
k = (kx^ky) is the wave vector parallel to the surface 
planes, and r(co’,k, d) is the total transmission probabil¬ 
ity of the electromagnetic waves. Notice that the second 
integral in Eq. (§ is carried out over all possible direc¬ 
tions of k and it includes the contribution of both prop¬ 
agating waves with k < uo/c and evanescent waves with 
k > ujjc^ where k is the magnitude of k and c is the 
velocity of light in vacuum. The transmission coefficient 
r(co’, k, d) can be expressed as [68] 

T{uj,'k,d)= (3) 

J Tr { [i - [i - , k<oj/c 

I iv{[7t2i k>oj/c ’ 

where q 2 = yJitP' jc? — is the 2 :-component of the wave 
vector in the vacuum gap and the 2x2 matrices are 
the reflections matrices characterizing the two interfaces. 
These matrices have the following generic structure 

j , (4) 

\ ' ij ij J 

where with a, /3 = s^p is the reflection amplitude for 
the scattering of an incoming a-polarized plane wave into 
an outgoing /3-polarized wave. Einally, the 2x2 matrix 
V appearing in Eq. ^ is defined as 

T) = [i-n2in23e^^^^^]-\ (5) 

Notice that this matrix describes the usual Eabry-Perot- 
like denominator resulting from the multiple scattering 
between the two interfaces. 

In Appendixes A and B we provide an alternative 
derivation of the central result of Eq. (§ that empha¬ 
sizes the non-reciprocal nature of our problem. More im¬ 
portantly, we show explicitly how the different reflection 
matrices appearing in Eq. (§ can be computed within 
a scattering-matrix approach for anisotropic multilayer 
systems. This approach provides, in turn, a natural 
framework to analyze different issues that will be cru¬ 
cial later on such as the nature of the electromagnetic 
modes responsible for the heat transfer. 

The result of Eqs. and ^ reduces to the well- 
known result for isotropic media first derived by Polder 
and Van Hove [7]. In that case, the reflections matrices 
of Eq. @ are diagonal and the non-vanishing elements 
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are given by 


S,s _ ^^ 

Qi + Qj 

ejQi + eiqj' 


(6) 

( 7 ) 


where qi = eiiJ^ jc? — is the transverse or 2 :- 
component of the wave vector in layer % and e^(co’) is 
the corresponding dielectric constant. Thus, the total 
transmission can be written as r(co’, k, d) = Tq{u:^ k, d) + 
rp(ci;, k, d), where the contributions of s- and p-polarized 
waves are given by 


Ta=s,p(^,k,d) = ( 8 ) 

/ (i-k2TP)(i-k2TP)/|7^c.P, k<uic 

I 4Im{r2i^}Im{r23^}e k>uj/c ’ 

where = 1 — Throughout this work 

we focus on the analysis of the radiative linear heat con¬ 
ductance per unit of area, d, which is referred to as the 
heat transfer coefficient. This coefficient is given by 


d(T, d) = lim 

AT^ 0 + 


Q(Ti =T +AT,T 3 = T,d) 
AT 


( 9 ) 


where T is the absolute temperature that we assume 
equal to 300 K throughout this work. Additionally, we 
define the spectral heat flux as the heat transfer coef¬ 
ficient per unit of frequency. In the following sections, 
we apply the general results presented here to different 
materials and magnetic field configurations. 


III. POLAR SEMICONDUCTORS: InSb 

The first obvious question to be answered is: In which 
materials can a magnetic field modify the NFRHT? Since 
the thermal radiation of an object is primarily deter¬ 
mined by its dielectric function, we need materials in 
which this function can be modified by an external mag¬ 
netic field, that is we need magneto-optical (MO) ma¬ 
terials. Focusing on room temperature experiments, the 
MO activity must be exhibited in the mid-infrared. Thus, 
doped semiconductors, where the MO activity is due to 
conduction electrons, are ideal candidates [69]. In these 
materials, one can play with the doping level to tune 
the plasma frequency to values comparable to the cy¬ 
clotron frequency at experimentally achievable magnetic 
fields, which is an important requirement to have sizable 
magnetic-induced effects in the NFRHT (see discussion 
below). Moreover, in semiconductors the NFRHT in the 
absence of field is typically dominated by surface electro¬ 
magnetic waves (both SPhPs and SPPs), which in turn 
are known to be strongly influenced by an external mag¬ 
netic field (69] [70]. Thus, it seems natural to expect a 
magnetic-field modulation of NFRHT in semiconductors. 

There is a variety of semiconductors that we could 
choose to illustrate our predictions. In this section we 


focus on InSb for several reasons. First, it is a po¬ 
lar semiconductor where the NFRHT in the absence of 
field is dominated by two different types of surface waves 
(SPhPs and SPPs), which allows us to study a very rich 
phenomenology. Second, InSb has a small effective mass, 
which enables to tune the cyclotron frequency to values 
comparable to those of the plasma frequency with moder¬ 
ate fields. Finally, InSb has been the most widely studied 
material in the context of magnetoplasmons and coupled 
magnetoplasmons-surface phonon polaritons. Thus, the 
magnetic field effect in the surface waves has been very 
well characterized experimentally inHii]. 


A. Perpendicular magnetic field: The realization of 
hyperbolic near-held thermal emitters 


Let us first discuss the radiative heat transfer between 
two identical plates made of n-doped InSb when the mag¬ 
netic field is perpendicular to the plate surfaces, i.e. 
H = i7,z, see Fig. [^a). In this case, the permittivity 
tensor of InSb adopts the following form m 


eiiH) -ie2{H) 0 
e{H) = I ie2iH) ei(if) 0 
0 0 £3 


( 10 ) 


where 

€i{H) = eoo 1 1 
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Here, eoo is the high-frequency dielectric constant, ujl 
is the longitudinal optical phonon frequency, ujt is the 
transverse optical phonon frequency, j (m*eoeoo) 

defines the plasma frequency of free carriers of density 
n and effective mass m*, F is the phonon damping con¬ 
stant, and 7 is the free-carrier damping constant. Finally, 
the magnetic field enters in these expressions via the cy¬ 
clotron frequency ujc = eHjm'". The important features 
of the previous expressions are: (i) the magnetic field in¬ 
duces an optical anisotropy (via the modification of the 
diagonal elements and the introduction of off-diagonal 
ones), (ii) there are two major contributions to the diag¬ 
onal components of the dielectric tensor: optical phonons 
and free carriers, and (iii) the MO activity is introduced 
via the free carriers, which illustrates the need to deal 
with doped semiconductors. In what follows we shall con¬ 
centrate in a particular case taken from Ref. [73], where 
eoo = 15.7, UJL = 3.62 x 10^^ rad/s, ujt = 3.39 x 10^^ 
rad/s, F = 5.65 x 10^^ rad/s, 7 = 3.39 x 10^^ rad/s, n = 
1.07 X 10^^ cm“^, m"jm = 0.022, and = 3.14 x 10^^ 
rad/s. As a reference, let us say that with these parame¬ 
ters uJc — 8.02 X 10^^ rad/s for a field of 1 T. Let us point 
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out that in this configuration, and due to the structure of 
the permittivity tensor, the transmission coefficient ap¬ 
pearing in Eq. (© only depends on the magnitude of the 
parallel wave vector, which considerably simplifies the 
calculation of the radiative heat transfer. 

Let us now briefly review the expectations for the heat 
transfer in the absence of magnetic field. As we show 
in Fig. [^b), the heat transfer coefficient features a large 
near-held enhancement for gaps below 1 fim. For d < 100 
nm this enhancement is largely dominated by p-polarized 
evanescent waves and the heat transfer coefficient in¬ 
creases as 1/(P as the gap decreases, which are two clear 
signatures of a situation where the heat transfer is dom¬ 
inated by surface electromagnetic waves. This can be 
further conhrmed with the analysis of the spectral heat 
hux, see Fig.j^c), which in the near-held regime is dom¬ 
inated by two narrow peaks that can be associated to 
SPPs (low-frequency peak) and SPhPs (high-frequency 
peak), as it will become evident below. Thus, the case of 
InSb constitutes an interesting example where two types 
of surface waves contribute signihcantly to the NFRHT. 
Let us now see how these results are modihed in the pres¬ 
ence of a magnetic held. 

In Fig. I^a) we show the heat transfer coefficient as a 
function of the gap size for diherent values of the perpen¬ 
dicular magnetic held. There are three salient features: 
(i) the far-held heat transfer is fairly independent of the 
magnetic held, (ii) in the near-held regime (below 300 
nm) the magnetic held suppresses the heat transfer by 
up to a factor of 3 (see inset), and (iii) by increasing 
the held, the heat transfer coefficient tends to saturate 
at around 6 T, although it is slightly reduced upon fur¬ 
ther increasing the held above 10 T (not shown here). 
The strong modihcation of heat transfer due to the mag¬ 
netic held is even more apparent in the spectral heat hux. 
As one can see in Fig. [^b), the magnetic held not only 
distorts and reduces the height of the peaks related to 
the surface waves, but it also generates a new peak that 
shifts to higher frequencies as the held increases. This 
additional peak appears at the cyclotron frequency and 
its presence illustrates the high tunability that can be 
achieved. Notice, for instance, that for a held of 6 T the 
thermal emission at the cyclotron frequency is increased 
by almost 3 orders of magnitude with respect to the zero- 
held case. 

To shed more light on these results it is convenient to 
examine the transmission of the p-polarized waves, which 
can be shown to dominate the heat transfer for any held. 
We present in Fig. this transmission as a function of 
the magnitude of the parallel wave vector, k, and the fre¬ 
quency for a gap d = 10 nm and diherent values of the 
magnetic held. As one can see, at low helds the trans¬ 
mission maxima are located around a restricted area of 
k and cj, clearly indicating that surface waves dominate 
the NFRHT. Notice also that their dispersion relation 
is modihed by the held, see Fig. ib). By increasing 
the held, those areas are progressively replaced by areas 
where the maximum transmission is reached for a broad 



FIG. 2. (a) Heat transfer coefficient for n-doped InSb as a 
function of the gap for diherent values of the magnetic held 
perpendicular to the plate surfaces. The inset shows the ra¬ 
tio between the zero-held coefficient and the coefficient for 
diherent values of the held in the near-held region, (b) The 
corresponding spectral heat hux as a function of the frequency 
(and wavelength) for a gap of d = 10 nm and diherent values 
of the perpendicular held. The solid lines correspond to the 
exact calculation and the circles to the uniaxial approxima¬ 
tion where the oh-diagonal terms of the permittivity tensor 
are assumed to be zero. 


range of /c-values and hnally, the surface waves are re¬ 
stricted to the reststrahlen band ujt < oj < ujl for the 
highest helds, see Fig. [^d). What is the nature of these 
magnetic-held-induced modes? 


To answer this question and explain all the results just 
described, it is important to realize that the oh-diagonal 
elements of the permittivity tensor do not play a major 
role in this conhguration. This is illustrated in Fig. |^b) 
where we show that the approximation consisting in set¬ 
ting 62 = 0 in Eq. (10) reproduces very accurately the 
exact results for the spectral heat hux for arbitrary mag¬ 
netic helds. This means that the polarization conversion 
is irrelevant and the plates ehectively behave as uniaxial 
media where their permittivity tensors are diagonal: e = 
diag[e^^,e^^,ezz], where = ei and = £ 3 . Within 
this approximation, which hereafter we refer to as uniax- 
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FIG. 3. The transmission coefficient for p-polarized waves as a function of the magnitude of the parallel wave vector and 
frequency for InSb and a gap of d = 10 nm. The different panels correspond to different values of the magnetic field that is 
perpendicular to the surfaces. The horizontal dashed lines separate the regions where transmission is dominated by surface 
waves (SPPs and SPhPs) or hyperbolic modes of type I and II (HMI and HMII). The white solid lines correspond to the 
analytical dispersion relation of the surface waves of Eq. (12). 


ial approximation, it is easy to compute the dispersion re¬ 
lation of the surface electromagnetic modes in our geome¬ 
try (see Appendix D). In the electrostatic limit /c ^ cj/c, 
the dispersion relation of these cavity modes is given by 


7 ^ 1 / I y ^XX ^ 

ksw = In ± yp 
n \ ^xx T y^xx/^ 


( 12 ) 


with the additional constraint that both exx and ezz must 
be negative. In the zero-field limit, this expression re¬ 
duces to the known result for cavity surface modes in 
isotropic materials [26]. As we show in Fig. see white 
solid lines, the dispersion relation of Eq. nicely repro¬ 
duces the structure of the transmission maxima in those 
frequency regions in which surfaces waves are allowed 
{^xxi^zz < 0). It is worth stressing that this dispersion 
relation describes in a unified manner both the SPPs that 
appear below the reststrahlen band and the SPhPs due 
to the optical phonons. More importantly, this disper¬ 
sion relation tells us that the magnetic field reduces the 
parallel wave vector of the surface waves and restricts the 
frequency region where they exist. Indeed, at high fields 
the SPPs disappear, while the SPhPs are restricted to 
the reststrahlen band. Fig. [^d). These two effects are 
actually the cause of the reduction of the NFRHT in the 
presence of a magnetic field. But what about the other 
modes that appear by increasing the field? Their nature 
can also be understood within the uniaxial approxima¬ 
tion. As we show in Appendix C, the allowed values for 
the transverse component of the wave vector inside th ese 
uniaxial materials are given by qo = \J^xx^"^ j for 
ordinary waves and q^ — fo^ 

traordinary waves. The dispersion of the extraordinary 
waves can be rewritten as 


-h kl 


T o 


,2 ’ 


(13) 


a dispersion that becomes hyperbolic when exx and ezz 
have opposite signs [48] . This is exactly what happens in 


our case in certain frequency regions at finite field. This 
illustrated in Fig. [^b-d), where we have indicated the 
hyperbolic regions defined by the condition exx^zz < 0- 
Notice that those regions correspond exactly to the areas 
where the transmission reaches its maximum for a broad 
range of /c-values. This fact shows unambiguously that 
our InSb plates effectively behave as hyperbolic materi¬ 
als. More importantly, and as it is evident from Fig. 
we can easily modify the hyperbolic regions by changing 
the field. Thus, we can change from situations where 
the hyperbolic modes (HMs) coexist with both types 
of surface waves to situations where the HMs dominate 
the NFRHT, which is what occurs at high fields, see 
Fig. id). Moreover, contrary to what happens in most 
hybrid hyperbolic metamaterials, we can have in a sin¬ 
gle material HMs of type I (HMI), where exx > 0 and 
ezz < 0, and HMs of type H (HMII), where exx < 0 and 
> 0, see Fig. ib-d). 

Let us recall that what makes HMs so special in the 
context of NFRHT is the fact that, as it is evident from 
their dispersion relation, they are evanescent in the vac¬ 
uum gap and propagating i nside the hyperbolic material 
for k > uj/c (HMI) or k > y^Je^\uj/c (HMII). Thus, they 
are a special kind of frustrated internal reflection modes 
that exhibit a very high transmission over a broad range 
of /c-values that correspond to evanescent waves in the 
vacuum gap [46]. As shown in Ref. [46], the number of 
HMs that contributes to the NFRHT is solely determined 
by the intrinsic cutoff in the transmission, which has the 
form r{uj^k) oc ex.p{—2kd) for k ^ uj/c. From this con¬ 
dition it follows that the heat flux due to HMs scales as 
l/(i^ for small gaps, as the contribution of surface waves. 
This explains why the appearance of HMs as the field 
increases does not modify the parametric dependence of 
the NFRHT with the gap size. Notice, however, that in 
spite of the high transmission of these HMs, their appear¬ 
ance does not enhance the NFRHT because they replace 
surface waves that possess even larger /c-values (notice 
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that the conditions of HMs and surface waves are mutu¬ 
ally excluding). Thus, we can conclude that the NFRHT 
reduction induced by the magnetic field is due to both 
the modification of the surface waves and their replace¬ 
ment by HMs that, in spite of their propagating nature 
inside the material, turn out to be less efficient transfer¬ 
ring the radiative heat in the near-held region than the 
surface waves. 

Let us point out that within the uniaxial approxi¬ 
mation, the heat transfer can be obtained in a semi- 
analytical form. In this case, the transmission coefficient 
is given by the isotropic result of Eq. ^ , where the re- 
hections coefficients adopt now the form 


r 


s,s 

21 


= r 


s,s _ 

23 ~ 


Q2 - go 

q2 + go 


r 


p,p 

21 


_ rrP’P 
— '23 


^xxQ 2 Qe 
^xxQ 2 Qe 


(14) 

(15) 


The uniaxial approximation is also useful to under¬ 
stand the high held behavior of the NFRHT. The ten¬ 
dency to saturate the thermal radiation as the held in¬ 
creases is due to the to the fact that the cyclotron fre¬ 
quency becomes larger than the plasma frequency and 
the last term in the expression of e^x = ^i, see Eq. & 
progressively becomes more irrelevant. Thus, the per¬ 
mittivity tensor becomes held-independent and the heat 
transfer is simply given by the result for uniaxial media, 
where Czz = has the form in Eq. & but = ei 
does not contain the last term in the hrst expression 
of Eq. (11). We hnd that the strict saturation of the 
NFRHT occurs at around 20 T and there is an interme¬ 
diate regime, between 6 and 20 T, in which the near-held 
thermal radiation slightly increases upon increasing the 
held (not shown here), leading to a nonmonotonic behav¬ 
ior. This behavior is due to an increase in the efficiency 
of the HMs that dominate the NFRHT in this high-held 
regime. 


To conclude this subsection, let us explain why the far- 
held heat transfer is rather insensitive to the magnetic 
held. For gaps much larger than the thermal wavelength 
(9.6 /im), the heat transfer is dominated by propagating 
waves and, as we show in Fig. the spectral heat hux 
in the absence of held exhibits a broad spectrum with a 
peak at around 1.5 x 10^^ rad/s. Indeed, the spectrum 
is very similar to that of a dielectric with a frequency- 
independent dielectric constant e = eool, see dotted line 
in Fig. As we illustrate in that hgure, the presence 
of a magnetic held only modihes this spectrum in a sig- 
nihcant way in a small region around the cyclotron fre¬ 
quency. This fact leads to a tiny modihcation of the heat 
transfer upon the application of an external held. As 
it can be seen in the inset of Fig. the magnetic held 
reduces the far-held heat transfer coefficient as the mag¬ 
netic held increases, but this reduction is quite modest 
and, for instance, it amounts to only 2.5% at a very high 
held of 20 T. 


Wavelength (pm) 



FIG. 4. Far-held spectral heat hux for InSb as a function 
of the frequency (and wavelength) for different values of the 
perpendicular held. These spectra have been computed for 
a gap d = 1 m. The dotted line corresponds to the result 
for plates made of a dielectric with a frequency-independent 
dielectric constant equal to Coo- The inset shows the corre¬ 
sponding ratio between the zero-held heat transfer coefficient 
and the coefficient for different values of the held. 


B. Parallel magnetic field 

Let us now turn to the case in which the magnetic 
held is parallel to the plate surfaces. For concreteness, 
we consider that the held is applied along the x-axis, 
H = iLa^x, but obviously the result is independent of 
the held direction as long as it points along the surface 
plane, as we have explicitly checked. In this case, the 
permittivity tensor of InSb adopts the form 

/ 63 0 0 \ 

e{H)= 0 -ie2 , (16) 

\ 0 ie 2 {H) ei{H) j 

where the e’s are given by Eq. (In}. Let us emphasize 
that in this case the transmission coefficient appearing 
in Eq. (§ depends both on the magnitude of the parallel 
wave and on its direction, which makes the calculations 
more demanding. Let us also say that we consider here 
the same parameter values for the n-doped InSb as in the 
example analyzed above. 

The results for the magnetic held dependence of the 
heat transfer coefficient for the parallel conhguration are 
summarized in Fig. [^a). As in the perpendicular case, 
the far-held is barely affected by the magnetic held, the 
near-held thermal radiation is suppressed by the held, 
and at high helds the NFRHT tends to saturates. Inter¬ 
estingly, it saturates to the same value as in the perpen¬ 
dicular conhguration. In spite of the similarities, there 
are also important differences. In this case, the NFRHT 
is much more sensitive to the held and a signihcant reduc¬ 
tion is already achieved at 1 T. Notice also that in this 
case the heat transfer coefficient is clearly nonmonotonic 













2 x10^ 10^ 2x10' 



FIG. 5. (a) Heat transfer coefficient for n-doped InSb as a 
function of the gap for different values of the magnetic field 
applied along the surfaces of the plates. The inset shows the 
ratio between the zero-field coefficient and the coefficient for 
different values of the field in the near-field region, (b) The 
corresponding spectral heat flux as a function of the frequency 
(and wavelength) for a gap oi d — 10 nm and different values 
of the parallel field. 


and the maximum reduction is reached at around 6 T. Fi¬ 
nally, notice also that the reduction is more pronounced 
than in the perpendicular case and the NFRHT can be 
diminished by up to a factor of 7 with respect to the zero- 
field case, see inset of Fig. |^a). This more pronounced 
reduction in the parallel configuration is also apparent in 
the spectral heat flux, as one can see in Fig. |^b). No¬ 
tice that also in this case there appears a high-frequency 
peak that is blue-shifted as the field increases. This peak 
appears at the cyclotron frequency and it has the same 
origin as in the perpendicular case. 

Again, to understand this complex phenomenology, 
it is convenient to examine the transmission of the p- 
polarized waves, which dominate the NFRHT for any 
field. Since in this case the transmission also depends 
on the direction of k, we choose to analyze the two most 
representative directions. In the first one, the in-plane 
wave vector k is parallel to the field, i.e. k = (/ca,,0), 
and in the second one, k is perpendicular to the field. 


i.e. k = (0,/c^). The transmission of p-polarized waves 
for these two directions is shown in Fig. |^as a function 
of the magnitude of the wave vector and as a function 
of the frequency for different values of the field. As one 
can see, the transmission exhibits very different behav¬ 
iors for these two directions. While for k || H the sit¬ 
uation resembles that of a perpendicular field (see dis¬ 
cussion above), for k T H it seems like the transmis¬ 
sion is dominated by surface waves that are severely 
affected by the magnetic field (with the appearance of 
gaps in their dispersion relations). These very different 
behaviors can be understood with an analysis of both 
the surface waves and the propagating waves inside the 
material in these two situations. In the case k || H, 
one can show that a uniaxial approximation, similar to 
that discussed above, accurately reproduces the results 
for the transmission found in the exact calculation. In 
this case, the permittivity tensor can be approximated 
by e = diag[exx,<^zz,<^zz\, where = es and = ei. 
Within this approximation, the dispersion relation of sur¬ 
face waves in the electrostatic limit /c ^ cj/c is also given 
by Eq. ([T^ (see Appendix D). As we show in Fig. HJa-c), 
this dispersion relation nicely describes the structure of 
the transmission maxima in the regions where the surface 
waves can exist (e^x, ^zz < 0)- On the other hand, as we 
show in Appendix C, the allowed values for the transverse 
component of the wave vecto r inside these uniaxial-like 
materials are given by Qo = for ordinary 

waves and — k‘^^xxl^zz foi* extraordi¬ 

nary waves. Again, the dispersion of these extraordinary 
waves is of hyperbolic type when €xx and €zz have oppo¬ 
site signs. In Fig.|^a-c) we identify the frequency regions 
where the HMs exist with the condition exx^zz < 0^ re¬ 
gions that progressively dominate the transmission as the 
field increases. Thus, we see that for k || H the situa¬ 
tion is very similar to that extensively discussed in the 
case in which the field is perpendicular to the materials’ 
surfaces. 

On the contrary, the situation is very different for 
k T H. In this case, there are no HMs and 
no uniaxial approximation can describe the situation. 
As we show in Appendi x C, the all owed g-values 
are given by go,i = and ^ 0,2 = 

(diy ^yz)^'^/~ which both describe waves 

with no hyperbolic dispersion. On the other hand, the 
dispersion relation of the surface waves in the electro¬ 
static limit is given by 

u — Jl ly. ~ ^ 'f'Vyz){Vyy — 1 — iVyz) \ /. 

2d [{pyy^l^iVyz){Vyy^l-iVyz)J^ ^ ^ 

where rjyy = tyy!{ely + and riy^ = -eyztiely + ejJ. 
As we show in Fig. I^d-f), this dispersion relation ex¬ 
plains the complex structure of the transmission maxima 
in this case. We emphasize that this dispersion relation 
is reciprocal in our symmetric geometry and for this rea¬ 
son we only show results for ky > 0. Notice that this 
dispersion is very sensitive to the magnetic field and al- 
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FIG. 6. The transmission coefficient for p-polarized waves as a function of the magnitude of the parallel wave vector and 
frequency for InSb and a gap of d = 10 nm. In all cases the field is parallel to the plate surfaces, H = The panels (a-c) 

correspond to different values of the magnetic field for wave vectors parallel to field, k = (kx,0), while panels (d-f) correspond 
to wave vectors perpendicular to the field, k = (0, ky). The horizontal dashed lines separate the regions where transmission is 
dominated by surface waves (SPPs and SPhPs) or hyperbolic modes of type I and II (HMI and HMII). The white solid lines 
correspond to the analytical dispersion relation of the surface waves of Eq. (12) in panels (a-c) and of Eq. 0 in panels (d-f). 


ready fields of the order of 1 T strongly affect the sur¬ 
face waves. Notice also the appearance of gaps in the 
dispersion relations, a subject that has been extensively 
discussed in the case of a single interface IMIIZO]- Over¬ 
all, the field rapidly reduces the /c-values of the surface 
waves and restricts the regions where they can exist. This 
strong sensitivity of the surface waves with k T H is the 
reason for the more pronounced reduction of the NFRHT 
for this field configuration. 

In general, for an arbitrary direction k = {k^^ky) the 
situation is somehow a combination of the two types 
of behaviors just described. The complex interplay of 
these behaviors for different k-directions is responsible 
for the nonmonotonic dependence with magnetic field, 
along with the change in efficiency of the HMs upon 
varying the field. On the other hand, at very high fields 
the cyclotron frequency becomes much larger than the 
plasma frequency and the off-diagonal elements of the 
permittivity tensor become negligible. At the same time, 
the field-dependent terms in the diagonal elements also 
become very small. Thus, the systems effectively become 
uniaxial and field-independent and the heat transfer is 
identical to the case in which the field is perpendicular. 
Finally, in the far-held regime, the heat transfer is not 
sensitive to the magnetic held for the same reason as in 
the perpendicular conhguration. 

Let us conclude this section with two brief comments. 
First, as it is obvious from the discussions above, an¬ 
other way to modulate the NFRHT is by rotating the 


magnetic held, while keeping hxed its magnitude. Actu¬ 
ally, we hnd that for any held magnitude, the NFRHT 
is always smaller in the parallel conhguration. Thus, one 
can increase or decrease the near-held thermal radiation 
by rotating appropriately the magnetic held. Second, we 
have focused here in the case of doped InSb, but similar 
results can in principle be obtained for other doped polar 
semiconductors such as GaAs, InAs, InP, PbTe, SiC, etc. 

IV. NON-POLAR SEMICONDUCTORS: Si 

In the previous section we have seen that when the held 
is parallel to the surfaces, one can have hyperbolic emit¬ 
ters, but the HMs always coexist to some degree with 
surface waves (even at the highest held). We show in 
this section that in the case of non-polar semiconduc¬ 
tors, where phonons do not play any role, it is possible 
to tune the system with a magnetic held to a situation 
where only HMs contribute to the NFRHT. For this pur¬ 
pose, we choose Si as the material for the two plates. 
As mentioned in the introduction, it has been predicted 
[431144] , and experimentally tested [23l |27] , that in doped 
Si the NFRHT in the absence of held can be dominated 
by SPPs even at room temperature. Let us see now how 
this is modihed upon applying a magnetic held. 

The dielectric properties of doped Si are similar to 
those of InSb, the only diherence being the absence of 
a phonon contribution. Thus, the dielectric functions of 
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FIG. 7. (a) Heat transfer coefficient for n-doped Si as a func¬ 
tion of the gap for different values of the magnetic field per¬ 
pendicular to the plate surfaces. The inset shows the ratio 
between the zero-held coefficient and the coefficient for differ¬ 
ent values of the held in the near-held region, (b) The cor¬ 
responding spectral heat hux as a function of the frequency 
(and wavelength) for a gap of d = 10 nm. 


Eq. now read 


^i(^) — ^oo ( 1 + 


Lol(uj + i'y) \ 

- (w + i7)2] j ’ 


£3 = Co 




(18) 


while € 2 {H) remains unchanged. Using the results of 
Ref. [43] for the dielectric constant of doped Si, we focus 
on a room temperature case where the electron concen¬ 
tration is n = 9.3 X 10^^ cm“^, Coo = 11-7, 7 = 8.04 x 10^^ 
rad/s, m*/m = 0.27, and ujp = 9.66 x 10^^ rad/s. We 
have chosen this doping level to have a situation in which 
the plasma frequency is not too high so that we can af¬ 
fect the NFRHT with a magnetic held, and not too low so 
that the NFRHT in the absence of held is still dominated 
by SPPs. 

The results for the heat transfer coefficient and spectral 
heat flux for a perpendicular magnetic field are displayed 


in Fig. Although there are several features that are 
similar to those of the InSb case, there are also some no¬ 
table differences. To begin with, notice that now higher 
fields are needed to see a significant reduction of the 
NFRHT (the required fields are around an order of mag¬ 
nitude higher than for InSb) and the reduction factors are 
clearly more modest, see inset of Fig. [^a). This is mainly 
a consequence of the smaller cyclotron frequency in the 
Si case for a given field due to its larger effective mass. 
Another consequence of the small cyclotron frequency is 
the fact that there is no sign of saturation of the NFRHT 
for reasonable magnetic fields. On the other hand, the 
spectral heat flux at low fields is dominated this time by 
a single broad peak that originates from SPPs (see dis¬ 
cussion below). As the field increases, the peak height is 
reduced and the peak itself is broadened and deformed. 
As we show in what follows, this behavior is due to the 
appearance of HMs that at high fields completely replace 
the surface waves. 


Again, we can gain a further insight into these results 
by analyzing the transmission of the p-polarized waves for 
different fields, which is illustrated in Fig. As one can 
see, the transmission is dominated by evanescent waves 
(in the vacuum gap) in a frequency region right below 
the plasma frequency. The origin of the structure of the 
transmission maxima can be understood with the uniax¬ 
ial approximation discussed above in the context of InSb. 
Again, this approximation reproduces very accurately all 
the results for arbitrary perpendicular fields (not shown 
here). Within this approximation, one can see that at 
low fields the transmission is dominated by SPPs, as we 
illustrate in Fig. [^a-b) in which we have introduced the 
dispersion relation of the SPPs given by Eq. ( 12 ). As 
soon as the magnetic magnetic field becomes finite, the 
system starts to develop HMs of type I in a tiny fre¬ 
quency region right above the region of existence of the 
SPPs, see Fig. |^b). The origin of these HMs is identi¬ 
cal to that of the InSb case, but the main difference in 
this case is that upon increasing the field, one reaches a 
critical field value (of 4.36 T for this example) for which 
the surface waves cease to exist and the transmission is 
completely dominated by HMs turning the Si plates into 
“pure” hyperbolic thermal emitters, see Fig. I^c-d). 


For completeness, we have also studied the heat trans¬ 
fer in the parallel configuration and the results for the 
heat transfer coefficient and spectral heat flux are shown 
in Fig. In this case the results are rather similar to 
those of the perpendicular configuration. In particular, 
contrary to the InSb case we do not find a nonmono¬ 
tonic behavior. Moreover, the NFRHT reduction is not 
much more pronounced than in the perpendicular case, 
although one can reach reduction factors of 50% for 12 
T. Finally, saturation is not reached for these high fields 
for the same reason as in the perpendicular configuration. 
As in the case of InSb, all these results can be understood 
in terms of the modes that govern the near-field thermal 
radiation. In this sense, for a direction where k || H, 
the SPPs that dominate the NFRHT at low fields are 
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FIG. 8. The transmission coefficient for p-polarized waves as a function of the magnitude of the parallel wave vector and 
frequency for Si and a gap of d = 10 nm. The different panels correspond to different values of the magnetic field that is 
perpendicular to the surfaces. The horizontal dashed lines separate the regions where transmission is dominated by surface 
plasmon polaritons (SPPs) or hyperbolic modes of type I (HMI). The white solid lines correspond to the analytical SPP 
dispersion relation of Eq. (p^ . 




CO (rad/s) 


FIG. 9. (a) Heat transfer coefficient for n-doped Si as a func¬ 
tion of the gap for different values of the magnetic field ap¬ 
plied along the surfaces of the plates. The inset shows the 
ratio between the zero-field coefficient and the coefficient for 
different values of the field in the near-held region, (b) The 
corresponding spectral heat hux as a function of the frequency 
(and wavelength) for a gap of d = 10 nm. 


progressively replaced by HMIs upon increasing the held 


and above 4.36 T they “eat out” all surface waves. On 
the contrary, for k T H there are no HMs and the only 
magnetic held effect is the modihcation of the SPP dis¬ 
persion relation. Again, the interplay between these two 
characteristic behaviors among the different k-directions 
explains the evolution of the NFRHT with the held. 

Let us conclude this section by saying that the behavior 
reported here for Si could also be observed for other non¬ 
polar semiconductors such as Ge. 


V. OUTLOOK AND CONCLUSIONS 

The results reported in this work raise numerous inter¬ 
esting questions. Thus for instance, in all cases analyzed 
so far, we have found that the magnetic held reduces the 
NFRHT as compared to the zero-held result. Is there 
any fundamental argument that forbids a magnetic-held- 
induced enhancement? In principle, there is no such an 
argument. The reduction that we have found in doped 
semiconductors is due to the fact that the we have ex¬ 
plored cases where surface waves, which are extremely 
efficient, dominate the NFRHT in the absence of held. 
In this sense, one may wonder if a held-induced enhance¬ 
ment could take place in a situation where the NFRHT in 
the absence of held is dominated by standard frustrated 
internal reflection modes, as it happens in metals m- 
Obviously, metals are out of the question due to their 
huge plasma frequency, but one can investigate non-polar 
semiconductors with a low doping level. Indeed, we have 
done it for the case of Si and again, we hnd that the mag¬ 
netic held reduces the NFRHT and moreover, exceed¬ 
ingly high helds are required to see any signihcant effect. 
Of course, we have by no means exhausted all possibili¬ 
ties and, for instance, we have not explored asymmetric 
situations with different materials. Thus, the question 
remains of whether the application of a magnetic held 
can under certain circumstances enhance the near-held 
thermal radiation. 
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The discovery in this work of the induction of hyper¬ 
bolic modes upon the application of a magnetic field may 
also have important consequences for layered structures 
involving thin films. Recently, it has been demonstrated 
that thin films made of polar dielectrics may support 
NFRHT enhancements comparable to those of bulk sam¬ 
ples due to the excitation of SPhPs [26]. Since hyperbolic 
modes have a propagating character inside the material, 
they may be severely affected in a thin film geometry by 
the presence of a substrate. Thus, one could expect much 
more dramatic magnetic-field effects in systems coated 
with semiconductor thin films. 

Obviously, the question remains of whether one can 
modulate the NFRHT with a magnetic field in other 
classes of materials. For instance, since a magneto¬ 
optical activity is required, what about ferromagnetic 
materials? Ideally, one could imagine to tune the 
NFRHT by playing around with the relative orienta¬ 
tion of the magnetization, following the spin-valve ex¬ 
periments in the context of spintronics. 

Another question of general interest for the field of 
metamaterials is if a doped semiconductor under a mag¬ 
netic field could exhibit the plethora of exotic optical 
properties reported in hybrid hyperbolic metamaterials 
|49l|50]. We have shown here that it can behave as a 
hyperbolic thermal emitter, but can it also exhibit neg¬ 
ative refraction or be used do to subwavelength imaging 
and focusing in the infrared? These are very important 
questions that we are currently pursuing. 

So in summary, we have presented in this work a very 
detailed theoretical analysis of the influence of a mag¬ 
netic field in the NFRHT. By considering the simple case 
of two parallel plates, we have demonstrated that for 
doped semiconductors the near-held thermal radiation 
can be strongly modihed by the application of an exter¬ 
nal magnetic held. In particular, we have shown that the 
magnetic held may signihcantly reduce the NFRHT and 
the reduction in polar semiconductors can be as large 
as 700% at room temperature. Moreover, we have shown 
that when the held is perpendicular to the parallel plates, 
doped semiconductors become ideal hyperbolic thermal 
emitters with highly tunable properties. This provides 
a unique opportunity to explore the physics of thermal 
radiation in this class of metamaterials without the need 
to resort to complex hybrid structures. Finally, all the 
predictions of this work are amenable to measurements 
with the present experimental techniques, and we are 
convinced that the multiple open questions that this work 
raises will motivate many new theoretical and experimen¬ 
tal studies of this subject. 
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Appendix A: Scattering matrix approach for 
anisotropic multilayer systems 

Our analysis of the radiative heat transfer in the pres¬ 
ence of a static magnetic field is based on the combina¬ 
tion of Rytov’s fluctuational electrodynamics (EE) and 
a scattering matrix formalism that describes the prop¬ 
agation of electromagnetic waves in multilayer systems 
made of optically anisotropic materials. As we show in 
Appendix B, the radiative heat transfer can be expressed 
in terms of the scattering matrix of our system. Thus, it 
is convenient to first discuss in this appendix the scatter¬ 
ing matrix approach employed in this work ignoring for 
the moment the fluctuating currents that generate the 
thermal radiation. Later in Appendix B, we show how 
this approach can be combined with EE. We follow here 
Ref. [76], which presents a generalization of the formal¬ 
ism introduced by Whittaker and Culshaw in Ref. HZ] 
for isotropic systems. 

Let us first describe the Maxwell’s equations to 
be solved. Assuming a harmonic time dependence 
exp(—icjt), the Maxwell’s equations for non-magnetic 
materials and in the absence of currents adopt the fol¬ 
lowing form: V • eocE = 0, V-H = 0, VxH = —icjeoeE, 
and V X E = icj/ioH, where the permittivity is in general 
a tensor given by Eq. 0 . The first Maxwell’s equation 
is automatically satisfied if the third one is fulfilled, and 
the second one can be satisfied by expanding the mag¬ 
netic field in terms of basis functions with zero diver¬ 
gence. Eollowing Ref. HZ], it is convenient to introduce 
the rescaling: cjeoE —> E and = ujjc ^ uj. Thus, 

the final two equations to be solved are 

V X H = -ieE, (Al) 

V X E = (A2) 

We consider here a planar multilayer system grown 

along the 2 :-direction in which the tensor e is constant 
inside every layer, i.e. it is independent of the in-plane 
coordinates r = (x, y). Thus, for an in-plane wave vector 
k = (kx^ky), we can write the fields as 

H(r, z) = h{zy'^'^ and E(r, z) = e{zy'^'^ . (A3) 

With this notation, Eqs. and ( |A 2 [ ) can be rewritten 
as 

ikyh^i^Z^ kyi^Z^ i ^ ^ €xj^j (^) ('^'^) 

3 

h'^{z) - ikxhz{z) = (^ 5 ) 

3 

ikxhyi^z^ '^kyhxi^z^ i ^ ^ (^), (''^^) 

3 
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and 


ikye^{z) - ey{z) = iuj‘^h^{z) (A7) 

e'^{z) - ikxCziz) = iuj‘^hy{z) (A 8 ) 

^kx^yi^z^ ikyCxi^z^ — iuj (^A9) 

where the primes stand for dz> 

Now our task is to solve the Maxwell’s equations for 
an unbounded layer. For this purpose, we write the mag¬ 
netic field h{z) as follows 


h{z) = |(/)a;X + 4>yy - ^{kx(j>X + ky<py)^ , (AlO) 

where x, y, and z are the Cartesian unit vectors and q 
is the 2 ;-component of the wave vector. Here, (j)x and 
(j)y are the expansion coefficients to be determined by 
substituting into Maxwell’s equations. Notice that this 
expression satisfies V • H = 0. Now, it is convenient to 
rewrite the previous expression in the vector notation: 


h{z) = ( 4>x,4>y,--{kx(t>x + ky(l)y) 


(All) 


With this notation, Eqs. ( A4| A 6 ) can be written as 


On the other hand, Eqs. ( A7||A9 ) adopt now the form 

e{z) = uj^h{z). (A13) 

Erom Eq. ( A12| ) we obtain the following expression for 
the electric field 


e{z) = fjCh{z)^ 


(A14) 


where 77 = e Substituting this expression in Eq. (A13) 
we obtain the following equation for the magnetic field 


C^fjCh{z) = uP‘h{z)^ 


(A15) 


which defines an eigenvalue problem for Indeed, only 
two of the three identities obtained from this equation, 
one for each x, y, and z, are independent. Erom the 
first two identities, and using Eq. (All), we obtain the 


following equations determining the allowed values for q 




(Alb) 


( 0 q —ky\ where (p = (px^^Py)^ and the 2 x 2 matrices An are defined 

—qO 1 . (A 12 ) by 
ky kx 0 J _I 


A 2 


y Vxy Vxx J 




kyVzZ 

^ -kxkyT^zZ 

_ I i^yi^x'lzx '^x^y'lzy '^y'lzx ^y^xyzy 

\ ^xVzy k^kyT]zx ^x^yVzy kykx'IJzx 


kxkyTjzz j I ^xVyy kxkyTJyx 
klv ZZ J \ kx kyTjxx ^xVxy 


1 0 
0 1 


(A17) 


This eigenvalue problem leads to the following quartic secular equation: = O 5 where the coefficients are 

given by 

D 4 = TjxxVyy Vxy'Hyx’) 

4^3 ~ kx \j]xyVyz T VyxVzy VyyiVxz T Vzx)] T ky IrjyxTjxz T VxyVzx VxxiVyz T Vzy)] : 

D 2 = k^ [qyy{Tlxx + Vzz) ~ VxyVyx ~ VyzVzy] + ky [TjxxiVyy T Vzz) ~ VxyVyx ~ VxzVzx] 

-\-kxky \jjxz{j]yz T Vzy^ T Vyz{Vzx Vxz)Vzz{Vxy T Vyx)] ^ {Vxx T Vyy)') 

= k^ [qxyVyz T VyxVzy VyyiVxz T Vzx)] T ky \VyxVxz T VxyVzx VxxiVyz T Vzy)] 

~^^x^y [VxyVzx T VxzVyx ~ VxxiVyz T Vzy)] T kykx [VyxVzy T VyzVxy ~ VyyiV XZ T Vzx)] 

-\-LU \k^{r]xz + Vzx) + ^y{Vyz + Vzy)] : 

Do = kl{r]yyr]zz - VyzVzy) + ky{r]xxVzz - VxzVzx) + k^ky [VxzVzy + VyzVzx “ Vzz{Vxy + Vyx)] 

~^kykx [VyzVzx T VxzVzy Vzz{Vyx T Vxy)] T k^ky [Vzzijlxx T Vyy) T VxyVyx VxzVzx VyzVzy] 

-\-UJ ^xiVyy Vzz) ky(^T]xx T Vzz) T kxky{7]xy T • (AI 8 ) 


In general, this secular equation has to be solved numeri- cally, but in many situations of interest the allowed values 
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for q can be obt ained analytically (see Appendix C). The 
solution of Eq. (A16) provides four complex eigenvalues 
for two lie in the upper half of the complex plane and 
the other two in the lower half. 

The next step toward the solution of the Maxwell’s 
equations in a multilayer structure is the determination 
of the fields in the different layers. This can be done 
by expressing the fields as a combination of forward and 
backward propagating waves with wave numbers q^ (with 
n = 1 , 2 ), and complex amplitudes and hn^ respec¬ 
tively. These amplitudes will be determined later by us¬ 
ing the boundary conditions at the interfaces and sur¬ 
faces of the multilayer structure. Since the boundary 
conditions are simply the continuity of the in-plane field 
components, we focus here on the analysis of the field 


components e^,, e^, and hy. From Eq. (|A 11 |), the 


in-plane components of h can be expanded in terms of 
propagating waves as follows 


hx{z) 

hy{z) 


t{( 

n=l ^ 


\ ^ (A19) 

J j 


where d is the thickness of the layer. Here, is the co¬ 
efficient of the forward going wave at the ^ = 0 interface, 
and bn is the backward going wave dX z = d. O n the other 
hand, correspond to the eigenvalues of Eq. (A16) with 
lm{qn} > 0 and Pn are the eigenvalues with Im{pn| < 0. 

To simplify the notation, we now define two 2x2 ma¬ 
trices and 4>_ whose columns are the vectors <pn and 
(Pn^ respectively. Moreover, we define the diagonal 2x2 
matrices f+( 2 ;) and t-{d — z), such that [t-^{z)]nn = 
and [f-{d — z)]nn = and the 2 -dimensional 

vectors h\\{z) = {hx{z),hy{z))'^ , a = (ai,a2)^, and 
b = ( 61 , 62 )^. In terms of these quantities, the in-plane 
magnetic-field components become 


hii{z) = 4 >+f+(^)a + 4 >_f_(d - z)b. 


(A 20 ) 


Similarly, from Eq. ( |A14 ) it is straightforward to show 
that the in-plane components of the electric field, ey (z) = 
{-ey{z),ex{z))'^, are given by 


611(2:) = ^ + .4^^$+ + A2^+q^ i+{z)a 

+ + A 2 ^-p) i-{d - z)b, (A 21 ) 


where the M’s are defined in Eq. (A17) and we have 
defined the 2 x 2 diagonal matrices q and p such that 

Qnn ~ Qn S^nd pnn ~ Pn- 

We can now combine Eq. (A20) and (A 21 ) into a single 
expression as follows 


V hl^) 


= M 


^ f+(2:)a \ 

f_((i — z)b J 


( Mil Mi2 \ ( f+(2)a 

VM21 M22 ) \L{d-z)b 


(A 22 ) 


where the 2 x 2 matrices Mij are defined as 

Mil = + Af^^+ + 

Mi2 = + A2^-P, 

M 21 — ^*+5 ^22 — 4*—. (A23) 

The final step in our calculation is to use the scatter¬ 
ing matrix (iS-matrix) to compute the field amplitudes 
needed to describe the different relevant physical quan¬ 
tities. By definition, the iS-matrix relates the vectors of 
the amplitudes of forward and backward going waves, ai 
and 6 /, where I now denotes the layer, in the different 
layers of the structure as follows 



Sii S 12 \ / ai' 

S 21 S 22 )\bi 


(A24) 

The field amplitudes in two consecutive layers are re¬ 
lated via the continuity of the in-plane components of 
the fields in every interface and surface. If we consider 
the interface between the layer I and the layer / + 1 , this 
continuity leads to 


f ^\\{di) 
V 


e||(0) \ 


(A25) 


where di is the thickness of layer /. From this condition, 
together with Eq. ( |A 22 ), it is easy to show that the am¬ 
plitudes in layers I and I + 1 are related by the interface 
matrix /(/, / -h 1 ) = in the following way 



= /(/,/ + !) 

= ( 

V ^21 I22 





(A26) 


where /,+ = ii,+ {di) and = f;+i _(di+i). 

Now, with the help of the interface matrices, the S- 
matrix can be calculated in an iterative way as follows. 
The matrix .§(/',/ + 1) can be calc ulate d from *§(/',/) 
using the definition of *§(/', 1) in Eq. (A24) and the inter¬ 
face matrix /(/,/ + 1). Eliminating ai and bi we obtain 
the relation between ai> ^ bi^ and a^+i, 6 /+ 1 , from which 
A§(/',/ + 1) can be constructed. This reasoning leads to 
the following iterative relations 


812 ( 1 ', I + 1) 


Ill-ftSl2{l',l)l2l] 


ill- 

X (ftSl2il',l)i22-il2)fr+l 


S2l{l', l + l) = S22{1', l)i2lSii{l', l + l)+ S2lil', 1) 
822 ( 1 ' ^ / + 1) = 822 ( 1 ' ^ 0^21^12(^^ / + 1) + 

^22(/',0^22/z+i. (A27) 


Starting from 8(1'J') = 1, one can apply the previous 
recursive relations to a layer at a time to build up 8(1' ,1). 
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Let us conclude this appendix by saying that from the 
knowledge of the ^'-matrix one can easily compute the 
field amplitudes in every layer and, in turn, the fields 
everywhere in the system ITT]. 


Appendix B: Thermal radiation in anisotropic 
multilayer systems 

In this appendix we show how the scattering matrix 
approach of Appendix A can be used to describe the ther¬ 
mal radiation between planar multilayer systems made of 
anisotropic materials. For this purpose, we first discuss 
how a generic emission problem can be formulated in the 
framework of the ^'-matrix formalism and then, we show 
how such a formulation can be used to describe the ther¬ 
mal emission of a multilayer system. 


1. Emission in the scattering matrix approach 


For concreteness, let us assume that there is a set of os¬ 
cillating point sources, with harmonic time dependence, 
occupying the whole plane defined by z = 2;'. The corre¬ 
sponding electric current density J is given by 

J(r, z) = 3o5{z - z') = joe^^'''5{z - z'), (Bl) 


where Jo(k) = Joe”**'’*'. This current density enters as 
a source term in Ampere’s law, Eq. ( |A1[), which now 
becomes V x H = J — ieE, while Eq. (A^ (Earaday’s 
law) remains unchanged. Thus, Eqs. (A4 adopt now 
the following form 

ikyhzi^Z^ ^ ^ ^ 

3 

h'^{z) - ik^h^iz) = joyS{z - z') -i'^eyjej{z) (B3) 

3 

ik^hy{z) - ikyh^{z) = joz^iz - z') - i'^e^jej{z). (B4) 


The presence of the source term induces discontinu¬ 
ities in the fields across the plane z = z' ^ as we proceed 
to show. Eirst, let us consider the effect of the in-plane 
components of the current density by putting jz = 0. To 
cancel the singular term due to the source in Eqs. (B2) 


and (B3), there must be discontinuities in and h 


dX z = z' equal to joy and —jox^ respectively. All the 
other field components are continuous, except for that 
exhibits a discontinuity equal to {kxjox + kyjoy)/€zz in 


virtue of Eq. (B4). Let us analyze now the role of the 
perpendic ular component of J by putting jox = joy = 0. 
Erom Eq. (B4), it is clear that in this case Cz must contain 
a singularity to cancel the singular term associated to 
the current source, that is ez{z) = —i{jQzl^zz)^{z — 
non-singular parts. This introduces sing ular term s in the 
left-hand side of the Maxwell Eqs. ( |A7[ ) and ( |A8[ ), which 
are cancelled by discontinuities in Cx and Cy equal to 
kxjoz/^zz and kyjo z/^z z^ res pect ively. Additionally, it is 
obvious from Eqs. ( |B2| ) and ( |B3[ ) that hx and hy acquired 
discontinuities equal to —^yz3{)zl^zz ^xzjoz!^zz^ re¬ 

spectively. Defining the following vectors 


P\\ — (jOy ^yzj^zl^zzi jox ^xzjoz/^zz) 

Pz = {-kyjozhzz,k:zjozl(^zzY', (B6) 


the boundary conditions on the in-plane components of 
the fields are thus 


)=Pz 

h\\{z'+) - = p\\. (B7) 

These discontinuity conditions can now be combined 
with the S'-matrix formalism of the previous appendix 
to calculate the emission throughout the system. Let 
us consider that the emission plane defines the interface 
between layers I and / + 1 in our multilayer structure. 
Thus, the boundary conditions in this interface become 


V ^11 (0) 



e\\{di) 



(B8) 


Using now the expression of thefields in terms of the 
layer matrices (M’s), see Eq. ( A22[ ), we can write 


M/+1 


t+i^i+i 




(B9) 


The external boundary conditions for an emission 
problem is that there should be only outgoing waves, that 
is ao = hw = 0, where 0 denotes here the first layer of 
the structure and N the last one. Using the defi nition s 
of the ^'-matrices 5(0, 1) and S{1 + 1, N) from Eq. ( |A24 ), 
it follows that 


ai = Si2{^,l)hi (BIO) 

= 52i(/ + l,7V)a,+i. (Bll) 


Substit utin g for ai and 6/+1 from Eqs. (BIO) and (Bll) 
in Eq. (B9) and rearranging things, we arrive at the fol¬ 
lowing central result 


[ Mii^/+1 + Mi2,Z+i/^^]^52i(^ + 1, A^) —[Mi 2 ,Z + Mii,//^^5i2(0,/)] \ / A _ f Pz 

\M2U+l+M22,l+lfr+lS2l{l + l,N) -[M22,; + M21.i/,+ 5i2(0,0] y V / VPlI 


(B12) 


which allows us to compute the field amplitudes on the left and on the right-hand side of the emitting plane. 
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From the solution of this matrix equation we can compute 
the field amplitude everywhere inside and outside the 
multilayer structure from the knowledge of the scattering 
matrix. 


2. Radiative heat transfer 

Let us now show that the previous results can be used 
to describe the radiative heat transfer. First of all, we 
need to specify the properties of the electric currents 
that generate the thermal radiation. In the framework 
of fluctuational electrodynamics [9], the thermal emis¬ 
sion is generated by random currents J inside the mate¬ 
rial. While the statistical average of these currents van¬ 
ishes, i.e. (J) = 0, their correlations are given by the 
fluctuation-dissipation theorem UEllTg] 

(Jfe(R,a;)J(*(R',w')) = ^^0(w,T)5(R - R')<5(a; - w') x 

TT 

[eM(R,w)-erfe(R,‘^)]/(20, (B13) 

where R = (r, and0(ci;,T) =/iC(;c/[exp(fico’c//cBF) —1], 

T being the absolute temperature. Let us remind the 
reader that with the rescaling introduced at the begin¬ 
ning of Appendix A, uj has dimensions of wave vector in 
our notation. Notice that in the expression of Q{uj^T) a 
term equal to hu:cj2 that accounts for vacuum fluctua¬ 
tions has been omitted since it does not affect the neat 
radiation heat flux. Notice also that we are using here 
the most general form of this theorem that is suitable 
for non-reciprocal systems. The fact that the current 
correlations are local in space and diagonal in frequency 
space reduces the problem of the thermal radiation to the 
description of the emission by point sources for a given 
frequency, parallel wave vector, and position inside the 
structure. Thus, we can directly apply the results derived 
in the previous subsection. 

Let us now consider our system of study, namely two 
parallel plates at temperatures T\ and T 3 separated by 
a vacuum gap of width d, see Fig. Our strategy to 
compute the net radiative heat transfer between the two 
plates follows closely that of the seminal work by Polder 
and Van Hove [7]. First, we compute the radiation power 
per unit of area transferred from the left plate to the 
right one, For this purpose, we first compute the 

statistical average of the ^-component of the Poynting 
vector describing the power emitted from a plane located 
at 2 ; = z' inside the left plate for a given frequency and 
parallel wave vector and then, we integrate integrate the 
result over all possible values of 2 :', cj, and k, i.e. 

pOO p pOO 

Qi^sid.Ti) = / duj dk dz'{Sz{oo,K z')). 

Jo J Jo 

(B14) 

A similar calculation for the power Qs^i transferred 
from the right plate to the left one completes the com¬ 
putation of the net transferred power per unit of area. 

Let us focus now on the analysis of the power emitted 
by a plane inside the left plate, see Fig.[^ This emitting 


0 

1 

2 

3 

^0 




Emitting^ 

plane 

z 

(left plate) 

d 

(vacuum gap) 

(right plate) 


FIG. 10. Two parallel plates separated by a vacuum gap of 
width d. The vertical dashed line inside the left plate indicates 
the position of an emitting plane that contains the radiation 
sources that generate the field amplitudes bo and ai. 


plane defines a fictitious interface between layers 0 and 
1, which are both inside the left plate. To determine the 
power emitted to the right plate we first compute the 
field amplitudes ai on the right hand si de of the plane. 
For this purpose we make of use of Eq. (B12), where in 


this case I = 0 and N = 3. Taking into account that 
'§ 12 ( 0 , 0 ) = 0 , it is straightforward to show that 

ai = — Mi2qM22|iM2iq] ^Pz + 

[M21,1 — ^22,1^12|i^11,i] ^P\\ 

= [Mi^iPz + [Mi^]i2P\\- (B15) 

To compute Qi^s, it is convenient to calculate the 
Poynting vector in the vacuum gap. For this purpo se, w e 
need the field amplitudes in that layer. From Eq. (A24), 
it is easy to deduce that these amplitudes are given in 
terms of ai as follows 


a 2 — T)§ii(l, 2 )ai, 

where D = [1 — § 12 ( 1 , 2 )§ 2 i( 2 , 3)]“^ and 

1 -1 


(B16) 


bo = 


1-521(2,3)^12(1,2) 


— § 21 ( 2 , 3 )a 2 . 


52 i( 2 , 3 ) 5 ii(l, 2 )ai 

(B17) 


It is worth stressing that the different elements of the 
scattering matrix that appear in the previous expressions 
can be factorized into scattering matrices S containing 
only information about the interfaces of the layered sys¬ 
tem, which are basically the Eresnel coefficients of the 
structure, and phase factors describing the propa gation 
between these interfaces. In particular, from Eq. (A27) 
it is easy to show that 


5ii(l,2) =5ii(l,2)/+(^') 

(B18) 

5i 2(1,2) =5i2(l,2)e*«^'^ 

(B19) 

52i( 2,3) =52i(2,3)e*«=^ 

(B20) 


where q 2 = is the 2 ;-component of the wave 

vector in the vacuum gap and 


/i+(^') = 


0 

o^q2,iz' 


(B21) 
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Here, Qi^i (with i = 1 , 2 ) are the 2 ;-components of the 
two allowed wave vectors in the medium 1 . On the other 
hand, the ^-matrices can be computed directly from the 
interface matrices as follows [see Eq. (A27)] 


5n(l,2)=7y(l,2) 

(B22) 

5i2(1,2) =-7^(1,2)112(1,2) 

(B23) 

^21 (2,3) =721(2,3)7^(2,3). 

(B24) 


In terms of the amplitudes a 2 and 62 , th e field s in the 
vacuum gap at z = 0 are given by [see Eq. (A 22 )] 


e||(0) 

ft.||(0) 


Mil ,2 [a2 - 

02 + e*«=‘'62 


(B25) 


where we have used that Mi 2,2 = —^ 11 , 2 , valid for any 
isotropic system. Thus, the ^-component of the Poynting 
vector evaluated at 2 ; = 0 in the vacuum gap reads 

S^{uj,k,z') = Ty^|/it(o)e||(0) + e|^(0)/i||(0)|^ 

- J_ /Zfo 

Auj Y eo 

2 _ a2 - 

(Mii,2 - ^> 2 } • (B26) 


|a2 (Mii,2 + -^11,2) ®2~ 
t 


Moreover, since 


Afii,2 — — 


1 f krK 


Q2 


\ =_ A 

kxky Uj‘^ — k‘^ 


1 

Q2 


(B27) 


and q 2 is either real (for k < uj) or purely imaginary (for 


k > uj)^ Eq. (B26) reduces to 

1 


5^(w,k, z') 


2q2UJ 
J 




(B28) 


a\Aa 2 — ^ 2 ^ 62 , 


k < UJ 


\ e '^^^^h\Aa 2 — e^^‘^^a\Ah 2 ^ k > uo, 

where the first term provides the contribution of propa¬ 
gating waves and the second one corresponds to the con¬ 
tribution of evanescent waves. 

Erom this point on, the rest of the calculation is pure 
algebra and we will not describe it here in detail. Let 


us simply say that the basic idea is to use Eqs. (B16) 
and (B17) to express the Poynting vector in Eq. ( B28D in 
terms of the field amplitude ai. Then, using Eq. (B15 ) 
and the fluctuation-dissipation theorem of Eq. (B13) one 


can calculate the statistical average of the Poynting vec¬ 
tor. Let us mention that the calculation can be greatly 
simplified by rotating every 2 x 2 matrix appearing in the 
problem from the Cartesian basis (x-y) to the basis of s- 
and p-polarized waves. This can be done via the unitary 
matrix 




kx ky 

ky kx 


(B29) 


which is the matrix that defines the transformation that 
diagonalizes the matrix A, i.e. 


Ad = RAR = 


0 q^2 


UJ" 0 

,2 


(B30) 


Einally, after in tegrat ing over all possible values of cj, k, 
and z', see Eq. (B14), one arrives at the following result 


for the power per unit of area transferred from the left 
plate to the right one 

Qi^3(rf,Ti) = ^~^e(a;,Ti) j ^T{co,k,d), 

(B31) 

where r(cc;, k, d) is the total transmission coefficient of the 
electromagnetic modes and it is given by 


r(co’, k, d) = 



[i — 5 'i2(1, 2)5'|2(1: 2)]D^[i — 5'2i(2, 3)5'2i(2, 3)]d| , k < uj (propagating waves) 

[>§ 12 ( 1 , 2 ) — Sl2{l^2)]D^[S2i{2^3) — § 21 ( 2 , 3)]T)| k > uj (evanescent waves) 


(B32) 


Here, the 2 x 2 matrices indicated by a bar are defined 
as follows 


D = 


(B33) 

(B34) 


Hollowing a similar reasoning, one can compute the 
power per unit of area transfer from the right plate to 


the left one and the final result reads 

Qz^i{d,TJ) = —e(w,T3) J j^^T{uj,k,d), 

(B35) 

where T{u},k,d) is also given by Eq. ( |B32 ). Thus, the 
net power per unit of area exchanged by the plates is 
given by Eqs. § and © in section [TT} To conclude, let 
us stress that in the manuscript uj is meant to be an 
angular frequency. 
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Appendix C: Dispersion relations 


amplitudes in the first and last layers as follows 


In this appendix we prov ide th e solu tion of the eigen¬ 
value problem of Eqs. (A16) and (A17) that give the dis¬ 


persion relations of the electromagnetic modes that can 
exist inside the materials considered in this work. In par¬ 
ticular, we focus here on three cases of special interest for 
our discussions in the main body of the manuscript. 

Case 1: e = diag[exxi^xxi^zz\- This situation is of 
relevance for the case in which the magnetic field is per¬ 
pendicular to the plate surfaces, see section |m|v. In this 
case, the allowed g-values are given by 


ql = 


ql = <^xx^‘^ - k'^ 




(Cl) 


Case 2: e = diag[exx^ ^zz]- This situation is rele¬ 
vant for the case in which the magnetic field is parallel 
to the plate surfaces, see section |IIIp , and the allowed 
g-values are given by 


Qo 


zOJ 


Qe 


c/e 


(C2) 


Case 3: the diagonal elements of e are exx and Cyy = 
ezzi while the only non-vanishing off-diagonal elements 
are eyz = —^zy This situation is relevant for the case in 
which the magnetic field is parallel to the plate surfaces, 
see section In this case, the allowed g-values adopt 
the following form 






^0,2 


_ / 2 I : 

— \eyy -h ey^ 




yy 


Appendix D: Surface electromagnetic modes 


ai 

bi 


■jv-i 




1=1 


dN 

bN 


^JS 


aN 

bN 


(D2) 

The condition for an eigenmode of the system is that 
ai = = 0 , which from the previous equation implies 

that ifiaN = 0. The condition for having a non-trivial 
solution of this equation is that det/f;^ = 0^ which is 
the condition that surface electromagnetic modes must 
satisfy. In our plate-plate geometry, the 4x4 matrix 
is simply given by 


/^ = J(l,2) 


e-i92di g 
6 


^(2,3), 


(D3) 


where let us recall that q 2 = Thus, the condi¬ 

tion for an eigenmode of the system reads 

det[Jii(l, 2)Jn(2, 2 )/ 2 i( 2 ,3)e*«='^] = 0. 

(D4) 

In what follows, we provide the explicit equations sat¬ 
isfied by the dispersion relation of the surface waves in 
the three cases considered in Appendix C. 

Case 1: In this case, Eq. (D4) leads to 


o-iq2d _ 


= ± 


Qe ^xxQ2 

Qe T ^xxQ2 


(D5) 


(C3) where is given in Eq. (Cl). This equations reduces to 
Eq. ([T^ in the electrostatic limit k ^ uj/c. 

Case 2: Here, assuming that the surface wave propa¬ 
gates along the x-direction, its dispersion relation satis¬ 
fies the following relation 


We briefly describe here how we determine the disper¬ 
sion relation of the surface electromagnetic modes in our 
system and we also provide the results for some configu¬ 
rations of special interest. 

Let us consid er a s tructure containing N planar lay¬ 
ers. Erom Eq. (A26), it is easy to show that the held 
amplitudes in layers I and I + I are related as follows 


= 


/t 0 


0 1 


-1 




1 0 

6 fi 


Z+l 


bi-\-i 




’/ + 1 


(DI) 


Now, using this relation recursively we can relate the held 


^-iq2d _ 


= ± 


Qe ^xxQ 2 
Qe T ^xxQ2 


(D6) 


where Qq is given in Eq. (C2). In the electrostatic limit, 
this equation reduces to Eq. (12). 

Case 3: In this case, and assuming that the surface 
waves propagate along the ^-direction, its dispersion re¬ 
lation is given by the solution of the following equation 


^-2iq2d _ 


{VyyQo,2 Q2 T 3yzk){gyyQo^2 Q2 Vyzk) 
{ilyyqo ,2 + q 2 + Vyzk){-qyyqo ,2 


where r]yy = eyy/{c 
go ,2 is given in Eq. 
equation reduces to 


q2 - Vyzk) ’ 

(D7) 
and 

C3|). In the electrostatic limit this 
;q. 


yy~l_^yzl Vyz ^yzj{_^ 


2 -fe2 , 

yy ^ ^yzJ 
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